home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / obsolete / betai.pro < prev    next >
Text File  |  1997-07-08  |  1KB  |  56 lines
  1. ; $Id: betai.pro,v 1.2 1997/01/15 04:02:19 ali Exp $
  2. ;
  3. ;  Copyright (c) 1991-1997, Research Systems Inc.  All rights
  4. ;  reserved. Unauthorized reproduction prohibited.
  5.  
  6. function betac, a,b,x
  7.  
  8.  lc = a + b
  9.  ln = a - 1.0
  10.  lq = a + 1.0
  11.  max = 100
  12.  ja = 1.0 & jn = 1.0 
  13.  
  14.  ka = 1.0 - lc*x/lq 
  15.  kn = 1.0
  16.  for i = 1,max DO BEGIN
  17.     ep = float(i)
  18.     tm = ep + ep
  19.     d = ep * (b-ep) * x/((ln + tm) * (a + tm))
  20.     jq = ja + d*jn
  21.     kq = ka + d * kn
  22.     d  = -(a + ep) * (lc + ep)*x/((lq + tm) * (a + tm))
  23.     jq1 = jq + d * ja
  24.     kq1 = kq + d*ka
  25.     prev = ja
  26.     jn = jq /kq1
  27.     kn = kq / kq1
  28.     ja = jq1 /kq1
  29.     ka = 1.0
  30.  
  31.     if ( ABS ( ja - prev) LT 3.0e-7 * ABS(ja)) then  return ,ja
  32.  
  33.  ENDFOR
  34.  
  35.  
  36. end
  37.  
  38.  function betai, x,a,b
  39.  
  40. if (x GT 1 OR x LT 0) THEN BEGIN
  41.   print, ' betat - x parameter out of bounds'
  42.   return, -1
  43.   ENDIF
  44.  
  45. ;gab =gamma(a) * gamma(b) 
  46. ;  gamma(a+b)/gab  * exp( a*alog(x) + b*alog(1.0-x))         $
  47.  
  48. if (x NE 0 and x NE 1 ) THEN temp =    $
  49.     exp(lngamma(a+b)-lngamma(a) -lngamma(b) + a*alog(x) + b*alog(1.0-x)) $
  50. else  temp = 0.0
  51.  
  52. if (x LT (a+1.0)/(a+b+2.0)) THEN return, temp * betac (a,b,x)/a     $
  53. else return, 1.0 -temp * betac(b,a,1.0 -x)/b
  54. END
  55.  
  56.